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Nonrelativistic electrostatic unmagnetized shocks are frequently observed in laboratory plasmas 
and they are likely to exist in astrophysical plasmas. Their maximum speed, expressed in units of 
the ion acoustic speed far upstream of the shock, depends only on the electron-to-ion temperature 
ratio if binary collisions are absent. The formation and evolution of such shocks is examined here 
for a wide range of shock speeds with particle-in-cell (PIC) simulations. The initial temperatures of 
the electrons and the 400 times heavier ions are equal. Shocks form on electron time scales at Mach 
numbers between 1.7 and 2.2. Shocks with Mach numbers up to 2.5 form after tens of inverse ion 
plasma frequencies. The density of the shock-reflected ion beam increases and the number of ions 
crossing the shock thus decreases with an increasing Mach number, causing a slower expansion of 
the downstream region in its rest frame. The interval occupied by this ion beam is on a positive 
potential relative to the far upstream. This potential pre- heats the electrons ahead of the shock even 
in the absence of beam instabilities and decouples the electron temperature in the foreshock ahead 
of the shock from the one in the far upstream plasma. The effective Mach number of the shock is 
reduced by this electron heating. This effect can potentially stabilize nonrelativistic electrostatic 
shocks moving as fast as supernova remnant (SNR) shocks. 

PACS numbers: 52.65.Rr, 52.35.Tc, 52.35. Qz, 98.38.Mz 



I. INTRODUCTION 



A supernova explosion accelerates a significant frac- 
tion of the material of the progenitor star up to a few 
percent of the hght speed c. The initial density of the 
blast shell plasma is many orders of magnitude larger 
than that of the surrounding plasma of the interstellar 
medium (ISM); the blast shell expands freely. The re- 
duction of the plasma density at the front of the radially 
expanding blast shell and its piling up of the ISM plasma 
imply that the latter eventually starts to affect the expan- 
sion. A forward shock between the upstream ISM plasma 
and the downstream plasma and a reverse shock between 
the downstream plasma and the (upstream) blast shell 
plasma form close to the blast shell front. The down- 
stream region between both shocks expands in time. The 
shock formation initiates the Sedov- Taylor phase of the 
supernova remnant (SNR) blast shell's expansion jj, in 
which it can propagate over astronomical distances be- 
fore it is stopped by its interaction with the ISM. 

The low particle collision frequency in the ISM implies 
that the forward shock is collision-less. Its dynamics 
is dominated by collective plasma processes, which im- 
poses constraints on the range of parameters for which 
stable shocks can form. The temperature of the cool 
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ionized component of the ISM is about 1 eV and the 
sound speed Cs « 10^ m/s. The sound speed is defined 
by Cg = (7sA:B/wi)(re -I- Ti) for electrons and ions with 
temperatures Tg and Ti (ion mass nii and Boltzmann 
constant ks), assuming that the adiabatic constant 7s is 
equal for both species. A fast SNR shock, like the north- 
eastern (NE) shock of the SNR RCW86 0, would have 
a Mach number Ms ~ 10^ - 10^ in the ISM plasma. The 
weak magnetization of the ISM 3] suggests that SNR 
shocks are essentially unmagnetized. 

Electrostatic unmagnetized plasma shocks can be gen- 
erated in laboratory plasmas 4, 5] and in simulation plas- 
mas [0, 0] by the collision of two plasma clouds at an 
appropriate speed. Collisions between identical plasma 
clouds can result in two types of unmagnetized electro- 
static shocks [3]- Sub-critical shocks can convert the 
entire kinetic energy of the inflowing upstream plasma 
into heat. The shock-reflected ion beam of super-critical 
shocks provides an additional energy dissipation mecha- 
nism and such shocks are stable at larger Mach numbers 
than the sub-critical shocks. Analytic estimates of the 
maximum Mach number exist, at which we find stable 
electrostatic shocks. The largest estimates of the maxi- 
mum Mach number of shocks, which develop out of the 
collision of two identical plasma clouds, are derived un- 
der the assumption that the electron temperature T^u 
exceeds by far the ion temperature Ti„ and that the par- 
ticle velocity distributions are Maxwellian's far upstream 
of the shock 0, Q • The peak Mach number of subcritical 



shocks is in this case Mg ~ 3, while that of super-critical 
shocks is Mg « 6.5. The peak Mach nuniber of shocks 
decreases with a decreasing ratio T^u/Tiu Q- No sta- 
ble electrostatic shocks exist under these approximations 
above a value Ms ~ 6.5, which is orders of magnitude 
below those of SNR shocks. 

An asymmetry between the colliding plasma clouds in 
terms of the electron temperature and density yields a 
double layer [9[ that can raise the maximum value of 
Ms [1(J]. Double layers are unstable, because the ther- 
mal pressure of the dense hot downstream plasma is not 
necessarily balanced by the deceleration or reflection of 
the upstream plasma. Experimental observations indi- 
cate that most double layers are rarefactive and some are 
compressive [ll| in their rest frame. They may thus not 
result in stable shocks, but they can trigger the forma- 
tion of stable electrostatic shocks [12] , even if the initial 
Mach number of the flow speed is 100 |13l |. 

Shocks can be stabilized for Ms > 6.5 by a magnetic 
field. An example is the Earth's quasi-perpendicular bow 
shock [1J|. It has a Mach number Ms ~ 10-12 with re- 
spect to the ion acoustic speed, if we take an electron tem- 
perature of 10 eV of the solar wind and a shock speed of 
500 km/s in the rest frame of the solar wind. This Mach 
number exceeds the stability limit of electrostatic shocks 
and the magnetic field must thus provide an additional 
stabilization. Indeed the Earth bow shock is a fast MHD 
mode shock |15| . Shocks that propagate into an upstream 
medium with a uniform magnetic field have been exam- 
ined in the context of SNR shocks with particle-in-cell 
(PIC) simulations [lfl-[2a|- These studies have addressed 
their stability, their efficiency in accelerating electrons 
and they have examined wave instabilities in the fore- 
shock [2a]. Usually a relatively strong magnetic field is 
used due to computational constraints. It is important 
to determine if a magnetic field is essential for a shock 
formation and for its stable propagation, given that the 
magnetic field of the ISM is relatively weak. 

Unmagnetized electrostatic shocks with a Mach num- 
ber Ms = 10^ may not exist in the ISM. This Mach 
number is however a least upper bound. It is based on 
the electron temperature of the ISM and neglects any 
pre-heating of the ISM by the shock. Unless the shock 
is quasi-perpendicular, shock-accelerated electrons or hot 
downstream electrons can escape upstream of the shock. 
They increase the ion acoustic speed, which in turn de- 
creases the Mach number. An electron foreshock is ob- 
served, for example, at the Earth's bow shock [23. Hot 
electrons are detected ahead of this shock even though 
the temperature difference between the electrons in the 
(downstream) magnetosheath and the upstream solar 
wind is well below that between the downstream region 
of SNR shocks and the ISM. Hotter electrons ahead of 
an SNR shock result in a higher Cg and therefore in a 
lower value of Ms . This pre-heating mechanism is inde- 
pendent of the heating of the foreshock by the cosmic ray 
precursor |28l - l30| and the wave precursors |31j . 

The possible existence of electrostatic unmagnetized 



SNR shocks motivates our numerical study. The shock 
is generated by letting an unmagnetized plasma collide 
head-on with a refiecting wall at the speed Vc- This cor- 
responds to a collision of two equal plasma clouds at 2vc, 
because the particles are reflected elastically. We per- 
form 7 separate simulations in which the collision speed 
between the plasma and the wall is varied from 1.06cs to 
2.48cs. The temperatures of the co-moving and spatially 
uniform electrons and ions are equal. The downstream 
plasma and its boundary, the shock, expand away from 
the wall and into the upstream region. The shock speed 
is thus always larger than the collision speed. 

Our results are as follows. All shocks are super-critical 
and move at a speed higher than 1.6 Cs in the upstream 
reference frame. The simulations thus demonstrate that 
stable shocks can form even if the ions are as hot as the 
electrons, which is not possible according to some ana- 
lytic models Q . A higher speed of the inflowing plasma 
results in a sharper shock discontinuity and in a larger 
density of the shock-reflected ion beam [7| . The fraction 
of ions that can cross the shock decreases and the expan- 
sion of the downstream plasma is slowed down; doubling 
Vc from 1.06 Cg to 2.12 c^ can only raise the shock's Mach 
number from 1.69 to 2.31. Shocks form instantly only for 
''^c ^ 1.9 Cs- An increase of Vc beyond this value delays 
the shock formation by tens of inverse ion plasma fre- 
quencies. A collision at Vc — 2.29cs results in a shock 
with Ms = 2.5 and no shock forms for Vc = 2.48cs dur- 
ing the resolved timescale, which is comparable to that 
achieved in laboratory experiments. It may thus not al- 
ways be possible to observe fast electrostatic shocks in 
laser-plasma experiments, while they can develop in the 
essentially unbounded astrophysical plasmas. 

The density of the shock-reflected ion beam, which in- 
creases with the collision speed, raises the overall ion 
number density ahead of the shock. This foreshock re- 
gion goes on a positive potential relative to the upstream 
plasma, which is maintained by the ambipolar electro- 
static held of the plasma density gradient. Upstream 
electrons are accelerated towards the shock by this elec- 
tric field. The accelerated electrons interact with the 
electrons, which leak from the downstream region into 
the foreshock. Their mixing raises the mean thermal en- 
ergy of the electrons in the foreshock to about 1/3 of the 
downstream one for all considered cases and the electron 
velocity distribution becomes a fiat-top one rather than a 
Maxwellian one. The simulations suggest that this elec- 
tron temperature depends on the relative speed between 
the pre-accelerated upstream electrons and the leaking 
downstream electrons and not on the electron tempera- 
ture far upstream of the shock. The implication is that 
once an electrostatic shock is present it is no longer ap- 
propriate to determine its stability properties using the 
Mach number it would have in the far upstream, because 
it propagates through a much hotter foreshock plasma. 

We outline the particle-in-cell simulation scheme and 
our initial conditions in section 2, we present our results 
in section 3 and we discuss them in section 4. 
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FIG. 1: The incoming plasma with density no moves to the 
left and collides with a reflecting boundary at a:: = 0. Initially 
the density of the reflected and incoming plasma adds up to 
2no. Electrons stream out of the overlap layer due to their 
larger thermal speed. A charge layer builds up at the front of 
the reflected plasma, which results in an electrostatic field. It 
confines the electrons and accelerates the ions to the right. 



We employ a ID simulation geom.etry, which separates 
the shock physics from, that of secondary instabilities and 
simplifies the interpretation of the results. Instabilities 
can be triggered by the counterstreaming ion beams in 
the foreshock with a relative speed that is larger than 
the ion acoustic speed. However, these instabilities can 
only develop if wave vectors are resolved that are ori- 
ented obliquely to the difference vector between the mean 
speeds of both ion beams [33l - t35j . Such vectors are ex- 
cluded here by the alignment of the collision velocity vec- 
tor with the simulation direction. The maximum shock 
speed is a few times the ion acoustic speed and thus 
much smaller than the electron thermal speed, which sup- 
presses the Buneman instability that can drive electro- 
static waves with a beam-aligned wave vector [36|. 



B. The solved equations and the initial conditions 



II. THE SIMULATION CODE AND THE 
INITIAL CONDITIONS 

A. The model 

The mechanism how a shock forms can be understood 
with the qualitative model shown in Fig. [TJ A plasma 
cloud is reflected at the position x=0. The cloud is spa- 
tially uniform and unmagnetized. It consists of electrons 
and ions, each with the number density no, that move at 
the speed modulus Vc towards x=0. We assume in this 
qualitative model that the ions are cold, so that the ion 
front does not spread out in time. The reflected ions and 
the incoming ions interpenetrate and a plasma sheet with 
an ion density 2no develops close to x = 0. It expands 
at the speed Vc to larger values of x. We call this sheet 
the overlap layer. The higher electron mobility implies 
that some electrons stream out of this layer. They form 
a thin negatively charged sheet just outside of it. An 
ambipolar electrostatic field is built up and the overlap 
layer goes onto a positive potential with respect to the 
inflowing upstream plasma. The electric fleld decelerates 
the ions that enter the overlap layer and this plasma com- 
pression raises the density of this layer above the value 
2no . This field also confines the electrons and accelerates 
the surface ions of the overlap layer [l^, IsJ . 

If the electric field is strong enough to equalize the 
speed of all ions in the overlap layer, then the ions of 
both plasma slabs will mix and form a single population. 
This downstream region is bounded by a forward shock 
and it expands due to ion accumulation. The upstream 
ions and electrons are heated up and compressed as they 
enter the downstream region. If all incoming ions enter 
the downstream region and mix with its plasma, then the 
shock is sub-critical. A fraction of the ions is refiected 
by a super-critical shock. They outrun the shock and 
form together with the incoming plasma the foreshock 
plasma. The upstream region is the one that has not yet 
been reached by shock-refiected ions. 



A particle-in-cell (PIC) simulation code [SJl solves 
Maxwell's equations together with the relativistic 
Lorentz force equation for an ensemble of computational 
particles (CPs). The collective charge distribution p{x) 
and current distribution J(x) on the grid is obtained from 
the interpolation of the charge and current of all CPs 
onto the grid. The collective charge and current distri- 
butions evolve in time the electromagnetic fields through 
the discretized forms of the Ampere's and Faraday's laws 
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They fulfiU Gauss' law V • E = p/eo and V • B ^ either 
as constraints or through correction steps. The electric E 
and magnetic B fields update in turn the momentum of 
each CP through the relativistic Lorentz force equation. 
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This equation updates the momentum p^ of the j 
tide of species i. The position i s up dated as dxj /dt 
We use the EPOCH PIC code ^. 

Our initial conditions are as follows. A plasma consist- 
ing of electrons and ions with mass ratio rrii/me = 400 
is introduced into a one-dimensional simulation box with 
length L. We use perfectly refiecting boundary condi- 
tions at x = and open boundary conditions at x = L. 
The spatially uniform plasma fills up the entire simu- 
lation box at t = 0. The particle number density of 
each species is uq and the initial temperature of their 
Maxwellian velocity distributions is Tq — 250ey. The 
electron thermal speed Ve = {kBTo/rrie) = 6.625 x 10^ 
m/s and Cs = 6.05 x 10^ m/s. The Debye length, 
which includes the contribution from the ions, is Xjj = 

Ve/{\/2ujp), where ojp = (noe^/eoWe) is the electron 
plasma frequency. The ion plasma frequency is tUpi = 



L0p/2Q. The box length L is resolved by 14400 cells of 
size A.X = Xd- We represent the electron species and the 
ion species by 2000 particles per cell, respectively. 

We perform 7 separate simulations with different col- 
lision speeds v^'. Vc — 1.06cs (Run 1), Vc = \A2ca (Run 
2), Wc = l.ncs (Run 3), Vc = 1.95cs (Run 4), Vc = 2.12cs 
(Run 5), Vc = 2.29c,. (Run 6) and Vc = 2.48cs (Run 7). 



III. SIMULATION RESULTS 

Non-relativistic electrostatic shocks are in their most 
basic form one-dimensional plasma structures. We as- 
sume that they are planar. Such shocks are fully de- 
scribed by the distribution of the electrostatic potential 
along the shock normal {x) direction and by the phase 
space distributions fe{x, Vx) and fi{x, Vx) of the electrons 
and ions, respectively. The shock's electrostatic poten- 
tial U{x, t) = — J^ E{x, t)dx has to be strong enough to 
reflect the incoming upstream ions, so we normalize it as 
2eU /rriiv'^. The potential difference between the far up- 
stream region and the downstream region should be > 1, 
because the shock speed exceeds Vc- Our reference po- 
tential is U{x = 0, t) at the reflecting boundary. We sub- 
tract from 2eU{x,t)/miv'^ the offset 2eUo/miv'^, which is 
the spatio-temporal average of the fully developed down- 
stream potential. The offset differs between the simula- 
tion runs, because we obtain in some cases electrostatic 
ion phase space holes [3^ at the wall, which put the wall 
and the downstream plasma on different potentials. 

The potential U — 2e{U {x ,t) — Uq) / rriivl reveals if and 
when the potential difference between the plasma over- 
lap layer and the incoming plasma cloud becomes large 
enough to reflect the incoming upstream ions. In what 
follows, we discuss separately the results of the simula- 
tion runs 1-7 followed by a comparison of the shocks in 
terms of plasma compression and electron heating. 

Run 1: The plasma collides at the speed Vc = l.OGcg 
with the wall and a shock should form, which is confirmed 
by the potential f/(a;, t) in Fig. [HJa). A positive potential 
forms at the wall on electron time scales. The front of 
this structure expands in the simulation frame at the 
constant speed v ~ 2000AD/(1200a;p/) or sa Q.&'ics for 
the ion-to-electron mass ratio that we use. The front 
of the potential thus expands at the speed ~ 1.69cs in 
the reference frame of the incoming (upstream) plasma. 
Typical values of the upstream potential reach « —1.8. 
The potential falls off smoothly in space and the width 
of the transition layer is of the order of lOOA^. 

The ion distribution in Fig. [DJb) shows that this tran- 
sition layer is actually wider. The plasma has reached a 
phase space distribution close to a; = 0, which is stable 
over the considered time scale. The velocity distribu- 
tion is not a single Maxwellian one that is centred at 
Vx — Q- We find instead a local minimum of the phase 
space density at u^; = and a; w with a value, which is 
about 20% below the peak one. The shock can thus not 
fully convert the flow energy of the incoming upstream 
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FIG. 2: (Color online) Collision speed Vc — 1.06cs (Run 
1): The distribution of the normalized electrostatic potential 
U{x,t) is shown in panel (a). The color scale is linear. The 
panels (b) and (c) show the phase space distributions of ions 
and electrons at the time tujpi — 600, respectively. The color 
scale is linear and normalized to the respective peak value. 



ions into heat. The width of the transition layer of the 
electrons in Fig. ^c) is even wider. It expands in time 
because the electrons react to the positive potential of 
the shock-reflected ion beam. The Buneman instability 
is suppressed by thermal effects and can not account for 
the observed electron temperature rise in the foreshock 
region. The electrons are heated instead by a turbu- 
lent mixing of electrons, which are accelerated by the 
foreshock potential towards the shock, with the electrons 
that leak from the downstream into the foreshock region. 
The downstream electrons (in the interval x < SOGAd at 
tupi — 600) do not have a Maxwellian distribution. The 
distribution fe{x « 0,^2:) shows a flat top centered at 
Vx = and a steep fall-off of the phase space density for 
ka;/'i^c| > 20. A flat-top distribution is also observed in 
the case of the foreshock electrons. 

The non-Maxwellian electron distribution alters the 
shock stability condition [4Q| and the ions do not get 
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FIG. 3: (Color online) Collision speed «c ~ 1.42cs (Run 
2): The distribution of the normalized electrostatic potential 
U{x,t) is shown in panel (a). The color scale is linear. The 
panels (b) and (c) show the phase space distributions of ions 
and electrons at the time tujpi = 600, respectively. The color 
scale is linear and normalized to the respective peak value. 



fully thermalized when they cross the shock. These two 
aspects may explain why a shock can form at the high ion 
temperature and Mach number we consider here, which 
exceeds significantly the limit determined by the analytic 
model in Rcf. [8]. 

Run 2: Raising the collision speed to Vc = 1.42cs 
leaves unchanged the qualitative structure of the shock. 
The potential distribution in Fig. [HJa) shows that 
the shock forms again on electron time scales and 
that it expands to larger x at a constant speed v w 
1600Ad/(1200Wp/), which corresponds to 0.52cs. This 
shock thus expands at a lower speed in the simulation 
frame than the one in Run 1. However, the larger Vc im- 
plies that its speed in the upstream frame of reference is 
higher with « 1.94cs. 

The ion distribution in Fig. ^h) shows some qualita- 
tive differences compared to the one in Fig. [DJb). The 
shock-reflected ion beam appears to be denser, which 
we confirm below. We observe in Fig. in^b) a pro- 



nounced velocity change of the incoming upstream ions 
at x/Xd ~ 800, which we do not find in this form in Fig. 
^b) . This stronger velocity change should lead according 
to the continuity equation to an increased plasma com- 
pression. The electron phase space distribution in Fig. 
^c) shows again a flat-top velocity distribution down- 
stream of the shock and a broad transition layer towards 
the upstream region. 

Figure [S^a) shows a weak depletion of the downstream 
potential that accelerates in time. This depletion crosses 
the position x/Xd ~ 350 at the time tujpi = 600 and it 
reaches x/Xu « 1300 at tojpi — 1200. This electrostatic 
field structure is connected to a depletion of the ion phase 
space density at x/Xd ~ 350 and v^/vc ~ +0.3 (Fig. 
[31[b)). We discuss ion phase space holes in more detail 
below. 

Run 3: The collision speed between the plasma and 
the wall is set to Vc/cs = 1.77. The maximum of the 
positive electrostatic potential in Fig. UJ^a) still develops 
close to the wall. The potential jump between the overlap 
layer and the upstream region does however not reach the 
value necessary for a shock formation on electron time 
scales as before. The potential jump reaches the value 
AC/ sa 1, which is required for a shock formation, at the 
time tujpi « 50. The front of the potential moves at the 
speed V w 600A_D/(600a;~/) or 0.39cs to larger x. Its 
Mach number in the upstream reference frame is 2.16. 
Figure mja) shows a potential structure that outruns the 
shock. This secondary structure reaches the position x — 
IOOOAd at tujpi ~ 200. Its speed is comparable to Vc, 
which corresponds to the speed at which the reflected 
ions move to larger x. The potential jump from U ~ 
to U sa —0.7 at the shock front after tujpi — 50, which 
separates the downstream region and the shock's foot 
(foreshock), is not enough to reflect ions that move at 
the speed —Vc towards the shock. However, the potential 
difference between the downstream region and the far 
upstream region is still large enough to reflect ions. 

Figure Hl^b) reveals that a shock is present at a; = 
600Ad. This structure reflects the incoming ions, which 
have the lowest speed in the reference frame of the shock, 
and it accelerates those ions from the downstream re- 
gion, which have the largest speed in the direction of the 
shock. The downstream population of the ions shows 
again a double-peaked velocity distribution that is get- 
ting more pronounced as we move away from the wall. 
This separation is upheld by the potential, which changes 
at tcdpi — 600 from the value C/ w at a; sa 200Xd to 
U ss —0.2 at a; ss 500 Ad- The shock transition layer 
has an extension comparable to a few hundred Xd- The 
electron distribution in Fig. \^c) shows again the subdi- 
vision into downstream electrons with x < 500Xd, fore- 
shock electrons in the interval 500 < x/Xd < 5000 and 
upstream electrons at even larger x. 

The potential distribution in Fig. IH^a) reveals a local- 
ized depletion of the potential in the downstream region 
that is similar to the one in Fig. Efa) but stronger. Its 
curved trajectory implies that it accelerates in time. It 
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FIG. 4: (Color online) Collision speed Vc — 1.77cs (Run 
3): The distribution of the normalized electrostatic potential 
U{x,t) is shown in panel (a). The color scale is linear. The 
panels (b) and (c) show the phase space distributions of ions 
and electrons at the time tujpi = 600, respectively. The color 
scale is linear and normalized to the respective peak value. 



reaches the shock position at the time tuipi = 600. This 
potential depression is caused by an ion phase space hole, 
which is demonstrated by Fig. [51 Ion phase space holes 
are structures, in which ions gyrate in a localized deple- 
tion of the electrostatic potential and form a phase space 
vortex. This negative potential is upheld by a local de- 
pletion of the ion charge density [33j ■ These vortices are 
tied to the ion distribution and move at speeds compara- 
ble to the ion thermal speed Vi . Its apparent acceleration 
in Fig. mja) arises from a change of the bulk speed of the 
ion beam that carries it. 

Run 4: A collision speed Vc/cg = 1.95 results in a 
shock formation that is similar to that in Run 3. The 
potential in Fig. [Hla) demonstrates that initially only a 
potential of C/ w —0.4 is reached close to the wall. A 
shock forms at tujpi sa 50, which propagates at the speed 
V w 500A£)/(600w~/) or w w 0.320^ to the right, yielding 
a Mach number Mg w 2.27 in the upstream frame of 
reference. A shock foot is visible that expands at a speed 



FIG. 5: (Color online) Panel (a) shows the electrostatic po- 
tential close to the shock at tupi = 300 (Run 3). Panel (b) 
displays the associated ion phase space density distribution 
on a linear scale. The potential depression at x/Xd ~ 200 in 
(a) is connected to the ion phase space hole in (b) at the same 
position and centered at v^/vc ~ 0.2. 



~ Vc to larger values of x. The ion distribution in Fig. 
^h) reveals a phase space vortex at a; ~ 0, a broad spatial 
interval 30 < a; < 500Ad with an ion distribution that is 
practically uniform along x and incoming and reflected 
ions at a; > 500. The electron distribution in Fig. [SJc) is 
qualitatively similar to those in the previous runs. 

Run 5: A simulation with Vc — 2.12cs shows a qualita- 
tively different shock formation stage compared to those 
in the previous cases. A weak potential U w —0.5 forms 
first. This structure detaches from the wall and its peak 
potential is located at x/Xd ~ 40 at tojpi « 70. This po- 
tential starts to grow and reaches [/ « at tcupi w 100. 
The potential decreases to C/ « —0.5 at tcUpi sa 120. We 
observe a cyclic reformation of the potential with a pe- 
riod of tujpi ~ 50. Despite this intermittent behaviour, 
the front of the potential propagates at a uniform speed 
400A£)/(600w~/) or 0.25cs into the upstream plasma. 
The front thus moves upstream at Ms « 2.37. The parti- 
cle distributions at tujpi = 600 in Figs. [7]Jb,c) agree with 
those observed in the previous two runs, implying that 
the shock structure does not strongly depend on the de- 
tails of its formation. An ion phase space hole has formed 
at the wall and the electrons have a flat top distribution 
downstream of the shock and in the foreshock. 

The cause of the initial oscillations of the potential is 
a non-stationary ion phase space hole at the reflecting 
boundary. Figure [8] compares the ion phase space dis- 
tributions close to the wall at the times tujpi = 100 and 



tw. 



pi — 120. The ion phase space hole and thus the mini- 



mum of the potential depression are centered at a; = at 
tujpi — 100. This implies that the plasma overlap layer 
is on a higher potential than the wall. The ion phase 
space hole's center and the minimum of the electrostatic 
potential have moved to a; ss IOXd at the later time. 
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FIG. 6: (Color online) Collision speed Vc — 1.95cs (Run 
4): The distribution of the normalized electrostatic potential 
U{x,t) is shown in panel (a). The color scale is linear. The 
panels (b) and (c) show the phase space distributions of ions 
and electrons at the time tojpi = 600, respectively. The color 
scale is linear and normalized to the respective peak value. 



FIG. 7: (Color online) Collision speed Vc — 2.12cs (Run 
5): The distribution of the normalized electrostatic potential 
U{x,t) is shown in panel (a). The color scale is linear. The 
panels (b) and (c) show the phase space distributions of ions 
and electrons at the time tujpi — 600, respectively. The color 
scale is linear and normalized to the respective peak value. 



The difference of the electrostatic potential between the 
plasma overlap layer and the wall has been reduced. 

Run 6: A collision speed Vc = 2.29cs is the largest one 
that results in a shock during a simulation time tupi = 
600. The shock does not form instantly, in line with 
the results of the simulation runs 3-5. It forms close 
to the wall and the potential distribution in Fig. M,^) 
shows no intermittent behaviour. The latter is thus not 
a consequence of large values oi Vc- A strong and steady 
foreshock potential is observed, which is again tied to 
the shock-reflected ion beam. We notice that a straight 
line fit of the shock front does not intersect x = at 
t = and the shock has thus not formed at the front 
of the plasma cloud. The ion distribution in Fig. ^h) 
agrees qualitatively with those in the simulation runs 3-5 
with respect to the subdivision into the ion phase space 
hole at the wall, a downstream region with a double- 
peaked velocity distribution and the foreshock structure 
consisting of incoming and reflected ion beams. 
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FIG. 8; Panel (a) and (b) show the ion phase space distri- 
bution at the times tujpi — 100 and 120, respectively. Both 
panels demonstrate that the ion phase space hole at the wall 
with a; < 30 is non-stationary. 
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FIG. 9: (Color online) Collision speed Vc — 2.29cs (Run 
6): The distribution of the normalized electrostatic potential 
U{x,t) is shown in panel (a). The color scale is linear. The 
panels (b) and (c) show the phase space distributions of ions 
and electrons at the time tojpi = 600, respectively. The color 
scale is linear and normalized to the respective peak value. 



An increase of the collision speed to Vc = 2A8Cs does 
not result in a shock formation during tujpi = 600, which 
is demonstrated by the corresponding U{x,t) displayed 
in Fig. [ini The potential difference between the plasma 
overlap layer and the far upstream region is not sufficient 
to yield a shock; no localized strong jump of the electro- 
static potential develops that could be associated with 
a shock. The ion phase space distribution (not shown) 
does not show a shock either. 

A comparison of the ion distributions in the simulation 
runs 1-6 shows two trends. The first trend is a shrinking 
of the ionic shock transition layer that connects the down- 
stream region with its spatially uniform double-peaked 
velocity distribution and the upstream region that con- 
sists of incoming and reflected ions. The width of this 
interval, which is characterized by maxima in the ion ve- 
locity distribution that diverge with increasing values of 
X, is 500Xd in Fig. [IJb) and it shrinks with increasing Vc 
to reach a few tens of Ad for the highest collision speed. 



FIG. 10: (Color online) The distribution of the normalized 
electrostatic potential U{x, t) computed from a simulation 
with Vc — 2.48cs is shown. The color scale is linear. 



The second trend, which has already been reported in 
Ref. [3], is that the density of the shock- reflected ion 
beam increases. The velocity distribution of the incom- 
ing ions and the reflected ions is practically symmetric 
with respect to fa; = for x > 300 Ad in Fig. ^h). 
The shock thus reflects practically all incoming ions and 
the low number density of ions that make it downstream 
explains the low expansion speed of 250Ad/(500w~/) or 
0.2cs in Fig. HJa). Its Mach number is M^ « 2.5. 

A comparison of the electron phase space distribu- 
tions obtained from the simulation runs 1-6 also reveals 
a trend. The distributions show a striking similarity if 
the velocity axis is scaled by Vc- An important feature is 
that the thermal spread of the electrons in the foreshock 
is practically unchanged by the choice of Vc- We attribute 
this to the way they are heated up. Figures mc)-[3I^c) and 
Figs. [Hc),[7jc),[9l^c) show that the electron's phase space 
distribution is bounded by two populations. The elec- 
trons with the largest positive values of Vc are supplied 
by downstream electrons that leak through the shock. 
We can determine the electron population that can leak 
into the foreshock as follows. The shock potential in the 
case of Run 6 can reflect ions that move from the up- 
stream region towards the shock at the speed Vc- The 
potential can thus confine downstream electrons with a 

1 /o 

velocity < [mi/m^) Vc- The fastest downstream elec- 
trons in Fig. Wic) have a speed modulus « 25fc and 
they can not be confined by the shock potential. They 
have lost a significant fraction of their kinetic energy by 
the time they have travelled far upstream but they are 
still fast in the foreshock region. The foreshock elec- 
trons with the largest negative speeds are supplied by 
the upstream electrons, which are accelerated towards 
the shock by the positive potential of the foreshock. Both 
counter-streaming electron populations thermalize in the 
foreshock. 



In what follows we want to provide further support for 
the trends outlined above. We analyse for this purpose 
the normalized ion density ni{x) = Uq J fi{x,Vx)dvx 
and the mean kinetic energy per electron Te{x) — 
neix)~ j vlfe{,x,Vx)dvx with ne{x) = J fe{x,Vx)dv^. 
We normalize T^{x) by the mean kinetic energy per elec- 
tron of a Maxwellian velocity distribution with the tem- 
perature Tq. 

We compare the ion density distributions in the Figs. 
ITlT a') for run 1-3 and ll^ a) for run 4-6. The correspond- 
ing mean thermal energies of the electrons are compared 
in the Figs. [TIT b) andfT^b). All curves are measured at 
the time tujpi — 600. We capture the full spatial interval 
of width QOOOXd affected by the shock-reflected ions. 

The ion densities downstream and, thus, the plasma 
compression increase with an increasing value of M^. 
This is a well-known fact that is here confirmed by the 
simulations. A somewhat surprising observation is that 
the downstream density in run 1 is only 2no, which cor- 
responds to a mere superposition of the density contri- 
butions of the incoming and reflected ions. The electro- 
static potential in Fig. ^a) was not strong enough to 
slow down significantly the incoming ions (See Fig. [2][b)) 
and the ions are thus not strongly compressed. 

The curves for rii (x) demonstrate that the width of the 
shock transition layer shrinks with an increasing value 
of Ms. The distribution ni{x) for run 1 changes from 
the downstream value to the foreshock value over several 
hundreds of Xd, while the ion density is practically dis- 
continuous on the displayed spatial scale in run 6. The 
ion density in the shock's foot increases with Mg. It is 
« 1.5no immediately ahead of the shock in run 1 and 
it gradually increases up to ~ 2no in run 6. The latter 
value implies that most ions are reflected by the shock; 
the shock acts as a piston Q. 

The curves for the mean thermal energy per electron 
follow the trend of the ion densities. The mean electron 
thermal energy is only slightly elevated in the foreshock 
region in run 1. The foreshock electrons get hotter with 
increasing Vc and the typical thermal energy of electrons 
in the foreshock has more than doubled in run 6 com- 
pared to the value far upstream. The shock in run 6 
thus moves through a plasma with an elevated sound 

1/2 

speed since Cs oc {Te +Ti) . This temperature rise is 
not large in our simulations. It can nevertheless be im- 
portant because it is caused by the mixing of leaking 
downstream electrons that move upstream and upstream 
electrons that have been accelerated towards the shock 
by the foreshock potential. The energy available to the 
incoming upstream electrons is dominated by the bulk ki- 
netic energy they gain in the foreshock potential and not 
by their thermal energy. The foreshock's electron tem- 
perature should thus be robust against changes of the 
electron temperature far upstream. The mean thermal 
energy of the electrons converges with increasing values 
of X to the initial value and the rise in the electron tem- 
perature is thus limited to the foreshock region. 

The mean thermal energies of the foreshock electrons 
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FIG. 11: (Color online) The total ion density in units of no 
is shown in panel (a). The mean thermal energy per electron 
normalized to that for To — 250 eV is shown in panel (b). 
The lowest (black) values corresponds to Vc/Cs — 1.06, the in- 
termediate (red) values to Vc/cs = 1.42 and the largest (blue) 



values to Vc/cs 



1.77. The time is tcjr, 



600. 
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FIG. 12: (Color online) The total ion density in units of no 
is shown in panel (a). The mean thermal energy per electron 
normalized to that for To = 250 eV is shown in panel (b). 
The lowest (black) values corresponds to Vc/cs — 1.92, the in- 
termediate (red) values to Vc/cs = 2.12 and the largest (blue) 



values to Vc/c 



2.29. The time is tui„ 



600. 



and the downstream electrons increase with Vc- We can 
compute from Figs. [TTT b) and fT^ b) the ratio of the mean 
electron thermal energies in the foreshock and down- 
stream. We compare for this purpose the respective val- 
ues at the positions x = ISOA^ and x = ISOOA^) in Fig. 
[TlTb) and those at a; = IOOAd and x = lOOOA^ in Fig. 
[T2jb). The measured ratios are displayed in Fig. [T3l as 
a function of the shock's Mach number. We find that 
the temperature of the foreshock electrons is in all cases 
about 1/3 of the value in the downstream region. 
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FIG. 13: The mean thermal energy of the foreshock electrons 
in units of the downstream value. 



IV. DISCUSSION 

We have performed here a series of one-dimensional 
particle-in-cell (PIC) simulations of electrostatic shocks. 
The purpose has been to gain insight into the conditions 
under which such shocks form and into the structure of 
their transition layer for a wide range of shock speeds. 
These aspects are important in the context of supernova 
remnant (SNR) shocks. Their large nonrelativistic speeds 
exceed the ion acoustic speed in the ISM by a factor 
10^ — 10^. This factor is much larger than the limit < 6.5 
for stable shocks, which is derived from various theoreti- 
cal models p-Q. It is thus essential to know the degree 
to which the ISM's temperature is important for SNR 
shocks. In other words, can a shock, once it is present, 
generate the electron temperatures it needs to be stable? 
This parametric study is also relevant for laboratory ex- 
periments. Important questions are what kind of shocks 
can form during the limited time over which observations 
are possible and what properties these shocks have. 

Our study addressed shocks in initially unmagnetized 
plasma that develop out of the collision of two equal 
plasma clouds at a speed that is less than 1% of the light 
speed c. Such shocks remain electrostatic because the 
low flow speeds yield only weak currents and magnetic 
fields. A speed in excess of 0.3c results in the growth 
of magnetic fields [2^, |3^, |41| that are strong enough to 
affect the shock dynamics and they become dominant at 
ultrarelativistic speeds even if the plasma was initially 
unmagnetized :42]. The chosen collision speeds are re- 
alistic for the late expansion phase of SNR blast shells. 
Our parametric study has provided the following results. 

Shocks form for the considered electron and ion tem- 
peratures up to a Mach number ^ 2.5. A shock with 
this Mach number reflects practically all incoming ions 
and acts as a piston. Raising the plasma collision speed 
by another 8% results in a shock-less interpenetration of 



the plasma clouds. Slower shocks with a Mach number 
Ms < 1.9 form on electron time scales while the faster 
ones develop slower, typically during about 10-100 in- 
verse ion plasma frequencies. The instability, which re- 
sults in the formation of shocks, does not set in immedi- 
ately after a collision of plasma clouds at a high speed and 
the formation time of the shock may be unpredictable. 
We can not exclude an eventual formation of a shock 
in our simulation that used the fastest collision speed. 
One possibility is that the ion acoustic instability, which 
can develop in multidimensional systems, pre-heats the 
plasma. The Mach number of the flow is reduced by 
the increasing ion acoustic speed and a stable shock may 
eventually form. Another shock formation mechanism is 
the destabilization of nonlinear plasma structures such 
as ion phase space holes that can trigger the formation 
of shocks even in very fast flows . 13| . 

The shock-reflected ion beam results in electron heat- 
ing in all simulations. This heating is not caused by 
instabilities such as the Buneman instability that is ex- 
cluded here by the low relative flow speed between elec- 
trons and ions. It arises from the electron acceleration in 
the ambipolar electric field of the plasma density gradi- 
ent. The foreshock electron temperature in our simula- 
tions has been about 1/3 of the electron temperature in 
the downstream plasma and it is probably independent of 
the upstream temperature. An electrostatic shock moves 
in this case through a medium with an ion sound speed 
that depends on the electron temperature in the post- 
shock plasma and on the density of the shock-reflected 
ion beam. The shock-reflected ion beam, which is essen- 
tial for the electron pre-heating, implies that the shock 
must be super-critical. 

The difference between the electron temperature in the 
foreshock and in the upstream is moderate in our simula- 
tions and at solar system shocks. The difference can be- 
come significant at astrophysical shocks, where the post- 
shock electron temperature can be of the order of keV 
while that in the upstream is about 1 eV. A foreshock 
temperature of SNR shocks that is a significant fraction 
of the downstream temperature implies that the Mach 
number of SNR shocks is in fact reduced by at least an 
order of magnitude compared to the least upper bound 
10^ — 10'^, which estimates the ion acoustic speed using 
the ISM temperature ~ 1 eV. 

The shock does not convert instantly the directed flow 
energy of the upstream into downstream heat. The shock 
potential and the ion distribution gradually change over 
a spatial interval that can be large if the Mach number is 
low, as our simulation shows. This has potentially impor- 
tant consequences for experimental observations of slow 
shocks. Their potential jump between the downstream 
region and the foreshock is moderate and it falls off over 
hundreds of Debye lengths. The electric field amplitude 
is thus low and smeared out over at least tens of Debye 
lengths. Slow shocks may thus not always be detectable 
in experiments that measure the electric field amplitude. 

Finally our simulations have shown that the ions down- 
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stream have a bi-Maxwellian velocity distribution and 
they are not in a thermal equilibrium. The shock- 
reflected ion beam is thus not the only process by which 
a supercritical shock can get rid of an excess flow en- 
ergy of the upstream plasma. This aspect together with 
the non-Maxwellian electron velocity distributions may 
help explaining why electrostatic shocks are stable even 
if electrons and ions have the same temperature, which 
is not possible according to some analytic models. 
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